function NumericalExperiments(model, parameter,filename,R,n,m,psi)

warning('off','all');

% memory preallocation
RPMLE = zeros(R,2); RSE_PMLE = RPMLE; RCI95_PMLE = RPMLE; RCI90_PMLE = RPMLE;
RGMM1 = zeros(R,2); RSE_GMM1 = RGMM1; RCI95_GMM1 = RGMM1; RCI90_GMM1 = RGMM1;
RGMM2 = zeros(R,2); RSE_GMM2 = RGMM2; RCI95_GMM2 = RGMM2; RCI90_GMM2 = RGMM2; 

r = 1;
while r<=R,
      % data generation
      [Y X] = DataGeneration(model,n,m,psi,parameter);
      % estimation
      [pmle se0 cond0] = PMLE(Y,X,psi); 
      [gmm1 se1 cond1] = GMMA(Y,X,psi); 
      [gmm2 se2 cond2] = GMMB(Y,X,psi); 
      % inference
      % confidence bounds at alpha = .05
      ci95_pmle = pmle-abs(norminv(.025))*se0 <= psi & psi <= pmle+abs(norminv(.025))*se0; 
      ci95_gmm1 = gmm1-abs(norminv(.025))*se1 <= psi & psi <= gmm1+abs(norminv(.025))*se1;
      ci95_gmm2 = gmm2-abs(norminv(.025))*se2 <= psi & psi <= gmm2+abs(norminv(.025))*se2;
      % confidence bounds at alpha = .10
      ci90_pmle = pmle-abs(norminv(.050))*se0 <= psi & psi <= pmle+abs(norminv(.050))*se0; 
      ci90_gmm1 = gmm1-abs(norminv(.050))*se1 <= psi & psi <= gmm1+abs(norminv(.050))*se1;
      ci90_gmm2 = gmm2-abs(norminv(.050))*se2 <= psi & psi <= gmm2+abs(norminv(.050))*se2;
      % writing away results
      convergence = cond0==0 & cond1==0 & cond2==0; 
      if convergence==1,
         RPMLE(r,:) = pmle'; RSE_PMLE(r,:) = se0'; RCI95_PMLE(r,:) = ci95_pmle'; RCI90_PMLE(r,:) = ci90_pmle';
         RGMM1(r,:) = gmm1'; RSE_GMM1(r,:) = se1'; RCI95_GMM1(r,:) = ci95_gmm1'; RCI90_GMM1(r,:) = ci90_gmm1';
         RGMM2(r,:) = gmm2'; RSE_GMM2(r,:) = se2'; RCI95_GMM2(r,:) = ci95_gmm2'; RCI90_GMM2(r,:) = ci90_gmm2';  
         
         r = r+1;
      end 
end

% saving output
save(filename);

median_r  = [   median(RPMLE);    median(RGMM1);    median(RGMM2)];
stddev_r  = [      std(RPMLE);       std(RGMM1);       std(RGMM2)]; 
ci95_r    = [mean(RCI95_PMLE); mean(RCI95_GMM1); mean(RCI95_GMM2)];
ci90_r    = [mean(RCI90_PMLE); mean(RCI90_GMM1); mean(RCI90_GMM2)];
display('--------------------------------------------------------------');
output = [median_r, stddev_r, ci90_r, ci95_r];   display(output);
display('--------------------------------------------------------------');

